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ABSTRACT 

In this letter we briefly describe the first results of our numerical study on the possi- 
bility of magnetic origin of relativist ic jets of long duration gamma ray bursters within 
the collapsar scenario. We track the collapse of massive rotating stars onto a rotating 
central black hole using axisymmetric general relativistic magnetohydrodynamic code 
that utilizes a realistic equation of state of stellar matter, takes into account the cool- 
ing associated with emission of neutrinos, and the energy losses due to dissociation of 
nuclei. The neutrino heating is not included. We describe the solution for one partic- 
ular model where the progenitor star has magnetic field B = 3 x lO^^G. The solution 
exhibits strong explosion driven by the Poynting-dominated jets whose power exceeds 
2 X 10^^ erg/s. The jets originate mainly from the black hole and they are powered via 
the Blandford-Znajek mechanism. 
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1 INTRODUCTION 

The phenomenon of Gamma Ray Burst (GRB) has been 
puzzling astrophysicists for many years since its discovery 
in 1970s (Klebesa del et al. 1973| |Muzets e t ari974 ). The 
recent identification of long duration GRBs with supernovae 
(see Delia Valle 2006, and Woosley & Bloom 2006 for full 
review) means that we are dealing with enormous amount 
of energy, 10^^ — lO^^erg, released within a very short time, 
2-100 seconds, in the form of highly relativistic collimated 
outflow ( Piran 2005 ) . Most of the current GRB studies are 



focused on the physics associated with production of gamma 
rays in such flows and their interaction with the interstel- 
lar medium or the stellar wind of the supernova progenitor. 
However, the central question in the problem of GRBs is 
undoubtedly the nature of their central engines. These pow- 
erful jets have to be produced as a result of stellar collapse, 
most likely by the relativistic object, neutron star or black 
hole (BH), formed in the center, and make their way through 
the massive star unscathed, remaining well collimated and 
highly relativistic. 

The most popular model of central engine is based on 
the "failed supernova" scenario of stellar collapse, or "col- 
lapsar", where the iron core of progenitor star forms a BH 
(Woosley 1993). If the progenitor is non-rotating then its 



collapse is likely to continue in a "silent" manner until the 
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whole star is swallowed by the BH. If, however, the specific 
angular momentum in the equatorial part of stellar enve- 
lope exceeds that of the last stable orbit of the BH then 
the collapse becomes highly anisotropic. While in the polar 
region it may proceed more or less uninhibited, for a while, 
the equatorial layers form dense and massive accretion disk. 
The gravitational energy released in the disk can be very 
large, more then sufficient to stop the collapse of outer layers 
and drive GRB outfiows, presumably in the polar direction 



where density is much lower ( ,MacFadyen &: Woosley 1999 ). 
In addition, there is plenty of rotational energy in the BH 
itself 



1/2^ 



(1) 



where Mbh is the BH mass and a G ( — 1, 1) is its dimension- 
less rotation parameter. For Mbh = BMq and a = 0.9 this 
gives the enormous value of E^ot — 8 x 10^^ erg. 

The three currently actively discussed mechanisms of 
powering GRB jets in the collapsar scenario are the heating 



via annihilation of neutrinos produced in the disk (Mac- 



Fadyen & Woosley 1999), the magnetic braking of the disk 
( Blandford Payne 1982| [Uzdensky MacFadyen 2006|), 
and the magnetic braking of the BH (Blandford & Znajek 



1977). The potential role of neutrino mechanism is rather 
difficult to assess as this requires accurate treatment of 
neutrino transport in a complex dynamic environment of 
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collapsar. The long and complicated history of numerical 
studies of neutrino-driven supernova explosions teaches us 
to be cautious. Numerical simulations by MacFadyen & 
Woosley( 1999) and Aloy et al.(2000) have demonstrated 
that sufficiently large energy deposition in the polar region 
above the disk may indeed result in fast collimated jets. 
However, the neutrino transport has not been implemented 
in these simulations and the energy deposition was based 
simply on expectations. When Nagataki et al.(2006) uti- 
lized a simple prescription for neutrino transport in their 
code they found that neutrino heating was insufficient to 
drive polar jets. 

A number of groups have studied the collapsar sce- 
nario using Newtonian MHD codes and implementing the 
Paczynski-Witta potential in order to approximate the grav- 



itational field of central BH (Proga et al. 2003; Fujimoto et 
al. 2006 Nagataki et al. 2QQ6| . In this approach it is impos- 
sible to capture the Blandford-Znajek effect and only the 
magnetic braking of the accretion disk can be investigated. 
The general conclusion of these studies is that the accre- 
tion disk can launch magnetically-driven jets provided the 
magnetic field in the progenitor core is sufficiently strong. 
Unfortunately, the jet power has not been given in most of 
these papers and is difficult to evaluate from the published 
numbers. In the simulations of Proga et al.('2003) the jet 
power at t 0.25s is 10^°erg/s. The initial magnetic field 
in these simulations is monopole with B c^i 2 x lO^^G at 
r — 3rg, where rg = GM^h/c^ (private communication). 

The study of collapsars in full GRMHD is still in its 
infancy. Sekiguchi & Shibataf2007') studied the collapse of 
rotating stellar cores and formation of BH in the collap- 
sar scenario. Their results show powerful explosions soon 
after the accretion disk is formed around the BH and the 
free falling plasma of polar regions collides with this disk. 
These explosions are driven by the heat generated as a re- 
sult of such collision. However, the authors have not ac- 
counted for the neutrino cooling and the energy losses due 
to photo-dissociation of atomic nuclei, and the explosions 
could be similar in nature to the "successful" prompt explo- 
sions of early supernova simulations (Bethe 1990). Mizuno 
et al.( |2004a[|2004b| carried out GRMHD simulations in the 
time-independent space-time of a central BH. The computa- 
tional domain did not include the BH ergosphere and thus 
they could not study the role of the Blandford-Znajek ef- 



fect (Komissarov 2004a). The energy losses have not been 
included and the equation of state (EOS) was a simple poly- 
trope. These simulations were run for a rather short time, 
280rg/c where rg = GM/c? ^ and jets were formed almost 
immediately due to unrealistically strong initial magnetic 
field. 

In this letter we describe the first results of axisymmet- 
ric GRMHD simulations of collapsars where we use realistic 
EOS ( Timmes k, Swesty 2000| , include the energy losses due 
to neutrino emission (assuming optically thin regime) and 
photo-dissociation of nuclei (see the details of micro-physics 
in Komissarov & Barkov 2007 ), use the computational do- 
main that includes the BH horizon and its ergosphere, and 
run simulations for a relatively long physical time, up to 
0.5s. The neutrino heating is not included. 



2 COMPUTER SIMULATIONS 

The simulations were carried out with 2D axisymmetric 
GRMHD code described in Komissarov ( 2004bt . Since this 
code can deal only with time-independent spacetimes we 
are forced to start from the point where the central BH has 
already been formed inside the collapsing star. In the pre- 
sented model the BH mass Mbh = BMq and its angular mo- 
mentum parameter a — 0.9. The mass density and the radial 
velocity of the collapsing star are described by the free-fall 
model of Bet he (1990 J corresponding to t = Is since the onset 
of collapse (see equations in Komissarov & Barkov, 2007). 
The parameter C is set to 9 corresponding to most massive 
stars. This gives us the free-fall mass rate M ^ O.5M0S~^. 
On top of this we endowed the free-falling plasma with an- 
gular momentum and poloidal magnetic field. The angular 
momentum distribution describes a solid body rotation up 
to the cylindrical radius — 6300 km. Further out the an- 
gular momentum is constant, / = lO^^cm^s"^. The magnetic 
field distribution is that of a uniformly magnetized sphere 
in vacuum, the radius of this sphere ri = 4500km and inside 
the sphere B — x 10^*^0. These distributions are intended 
to describe the progenitor at the onset of collapse rather then 
at the state developed one second later. We utilize the Kerr- 
Schild coordinates of spacetime. The computational grid is 
uniform in polar angle, ^, where it has 180 cells. In the ra- 
dial direction it is uniform in logr, and has 450 cells. The 
inner boundary is located just inside the event horizon and 
adopts the free-flow boundary conditions. The outer bound- 
ary is located at r = 25000/cm and at this boundary the flow 
is prescribed according to the Bethe's model. 

At the beginning of simulations the angular momentum 
of accreting gas is less than that of the last stable orbit, /iso- 
It falls straight into the BH, and no disk is formed. How- 
ever, the magnetic flux threading the BH gradually increases 
and so is the magnetic pressure. When the outer layers with 
/ > /iso reach the BH the centrifugal force slows down their 
infall and the accretion disk is beginning to form around the 
BH. At the same time the accretion shock separates from its 
surface. The low angular momentum plasma of polar regions 
k.png falling straight into the BH after passing the accretion 
shock whereas the high angular momentum plasma fills the 
"bubble" above and below the disk (fig[l]). Strong differen- 
tial rotation within this bubble leads to amplification of the 
azimuthal component of magnetic field, the magnetic pres- 
sure grows and eventually overwhelms the ram pressure of 
free- falling envelope - the explosion begins. The BH is a key 
player in the process pumping electromagnetic energy into 
the bubble and the disk at the rate of ~ 2 x lO^^ergs"^. 

Figures [2] and [3] show the solution at t = 0.45s, near the 
end of simulations. At this time, the solution exhibits two 
well defined polar jets surrounded by the magnetic cocoons 
of high pressure and low density. The magnetic pressure of 
these cocoons, which have been inflated by the jets, exceeds 
by more than six orders of magnitude the magnetic pres- 
sure in the collapsing star. These over-pressured cocoons 
drive strong bow shock (blast wave) into the star (right 
panel of fig[2]). The mean propagation speed of the shock 
in the polar direction ^ 0.18c. Near the equator the 
stellar plasma compressed by the shock continues stream- 
ing downward with supersonic speed. At the equator and 
well outside of the accretion disk the stream coming from 
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Figure 1. Solution immediately before the explosion (t=0.24s). Left panel: the baryonic rest mass density, logiop, in gjcvn? and the 
magnetic field lines; Middle panel: the ratio of gas and magnetic pressures, logioP/ Pm, and velocity direction vectors; Right panel: the 
ratio of azimuthal and poloidal magnetic field strengths, logioB'^ / Bp, and the magnetic field lines. 




Figure 2. Solution on different scales at f = 0.45s. The colour images show the baryonic rest mass density, logiop in g/cm^, the contours 
show the magnetic field lines, and the arrows show the velocity field. 



northern hemisphere collides with the stream coming from 
the southern hemisphere and a pair of oblique shocks de- 
velop at r ~ 200rg (middle panel of fig|3|. These shocks are 
not strong enough to cause photo-dissociation of nuclei and 
the high post-shock pressure drives the reflected flows away 
from the equatorial plane. Plasma from the skin layers of the 
reflected streams actually enters the bubbles and interacts 
with the jets (We expect this effect to weaken later when 
the blast wave moves further away.) The inner layers of the 
reflected streams pass through another shock, at r ^ 50r^, 
and feed the accretion disk. The left panel of fig[2] shows 
the solution in the immediate vicinity of the BH. Its struc- 
ture is reminiscent to that found in the previous studies of 
thick disks around BHs - main disk, its dynamic corona, 
and magnetically-dominated fun nel (|De Villiers Hawley| 



2003||McKinney k Gammie 2004| [Shibata et al. 2007| ). This 



funnel is the region there the Pointing dominated jets are 
produced as well as the "wind" blowing into the BH - in 
this image one can clearly see the surface separating these 
flows. 

Figure |3] shows the magnetic properties of the cen- 
tral region. Not only the funnel but also the disk corona 
are magnetically-dominated. The magnetic field strength 
reaches fewxlO^^G near the BH, it is weaker in the funnel 
compared to the disk at the same spherical radius but not by 
much. Within the disk and corona the azimuthal magnetic 
field, B"^ exceeds the poloidal one. Bp, by two or three orders 
of magnitude. In contrast, in the funnel B^ / Bp < 1, reaching 
unity only near the funnel walls. In fact, the poloidal field in 
the funnel exceeds that in the disk and corona by 1-2 orders 
of magnitude. This is in contrast to the conclusion made by 



Ghosh & Abramowicz( 1997) and Livio et al.(1999), namely 
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Figure 3. The inner region t = 0.45s. Left panel: the magnetization parameter, logioP/ Pm, and the magnetic field lines; Middle 
panel: the ratio of azimuthal and poloidal magnetic field strengths, logioB^ / Bp, and the magnetic field lines; Right panel: the magnetic 
field strength, logio(B), and the magnetic field lines. 



that the poloidal field threading the BH horizon should be of 
the same order as the poloidal field in the inner parts of the 
disk. Their main argument, that both fields are produced by 
the same azimuthal current flowing in the disk, misses the 
fact that additional currents may flow in the magnetosphere 
and over the disk/funnel surface and support the magnetic 
field inside the funnel in the manner similar to solenoid. In 
our case, the poloidal field threading the BH is the original 
field of the progenitor that has been accumulated during the 
initial phase of free infall. 

The left panel of fig|4] shows the baryonic mass flux 
as a function of spherical radius. One can see that it re- 
duces from the free-fall value M —0.5Mqs~^ down to 
M — O.OGMqs"^ at the event horizon. Between r 60r^ 
and r = 2500r^ this reduction reflects the effect of the bow 
shock driven into the star by the jets. The sharp reduction at 
r ^ 60rg corresponds to the position of the accretion shock 
and marks the transition from approximate free-fall to the 
centrifugally supported disk. 

The middle panel of fig|4] shows the integral energy 
fluxes of the jets as functions of spherical radius. To be 
more precise the integration is carried out over the the whole 
sphere but the contribution from areas with the baryonic 
rest mass density p > lO^gcm"^ is excluded. We have ver- 
ified that the bulk contribution to the fluxes computed in 
this way comes from the jets. The baryonic rest mass flux, 
pu^ radial component of 4- velocity, is excluded from the to- 
tal and the matter energy fluxes, that is these fluxes are 
computed via 



E = -27T 



Js 



+ pu)^de, 



(2) 



where 7 is the determinant of the metric tensor of space 
and T is either the total stress-energy-momentum tensor 
or its hydrodynamic part. The most important conclusion 
suggested by this figure is that at least 80% of the jet en- 
ergy is provided directly by the BH and at a very high rate, 
£^ ~ 2 X lO^^ergs"^. The remaining 20% seem to be pro- 
vided by the inner part of the disk - this explains the rise of 



jet power between the event horizon and r lOr^. Indeed, 
careful examination of the solution shows that some mag- 
netic field lines enter the jet from the skin layers of the disk 
with p > lO^gcm"^. However, it remains to be shown that 
this is not caused by the numerical diffusion of magnetic flux 
from the funnel into the disk. The right panel of fig^ shows 
the distributions of Poynting flux and hydrodynamic energy 
flux (including the rest mass-energy) across the horizon and 
allows us to determine whether it is the Blandford-Znajek 



or the MHD-Penrose mechanism (Punsly & Coroniti 1990 
iPunsly 2001 1 |Koide et al.2002t or both of them that pro- 



vide the energy supply to the jets. Since the hydrodynamic 
flux is everywhere negative the MHD-Penrose mechanism 
can be ruled out with certainty. This is confirmed by the 
fact that the hydrodynamic energy-at-infinity is positive ev- 
erywhere inside the ergosphere. Thus the jet is powered by 
the Blandford-Znajek mechanism. For a force-free monopole 
magnetosphere the Blandford-Znajek power is given by 



1 ([^Y 

6c V 47r y 



6c V 47r ' ' 

where Qh is the angular velocity of the BH and ^ is the mag- 
netic flux threading the BH. In the derivation we assumed 
that the angular velocity of magnetosphere Q = 0.5^^. This 
holds well even for rapidly rotating BHs with monopole mag- 
netospheres (Komissarov 2001 ) and corresponds to the mean 
value of Q measured in our simulations as well. Using the 
measured value of ^ we derive E^z — 2.6 x lO^^ergs"^ which 
agrees quite well with the value of Ebz provided by figji] The 
total amount of free energy-at-infinity in the bow shock and 
the bubble at time t = 0.45s is ^ 1.37 x lO^^erg. Since 
the explosion develops only at t = 0.24s the mean jet power 
over the active period is < ^ 6 x lO^^ergs"^, indicating 
the higher jet power at the early stages of the explosion. 

The middle panel of fig|4] also shows that initially the 
jets are Poynting-dominated but gradually the electromag- 
netic energy is converted into the energy of matter. However, 
the accuracy of our simulations is insufficient to capture the 
jet dynamics. First of all, we are forced to keep the flow 
magnetization below the limit at which the code crashes - 
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Figure 4. Left panel: the integral baryonic mass flux in units Mqs"-*^ as a function of spherical radius; Middle panel: the integral fluxes 
of total energy (solid line), electromagnetic energy (dashed line), and hydrodynamic energy (dotted line); Right panel: the energy flux 
densities at the event horizon for electromagnetic energy (solid line) and hydrodynamic (matter) energy. Time t = 0.45s. 



this is done via artificial injection of plasma in the danger 
cells. This reduces the length scale for the energy conversion 
via magnetic acceleration of plasma, as well as the asymp- 
totic Lorentz factor. The numerical mass diffusion into the 
jets from the disk corona further exacerbates this problem. 
Finally, numerical resistivity causes dissipation of the jet 
electric current. Due to the mass diffusion and numerical vis- 
cosity the jets never become ultrarelativistic - their Lorentz 
factor rarely exceeds F = 3. On the other hand, the total 
energy is conserved we do not expect these numerical prob- 
lems to have strong effect on the dynamics of the bow shock 
and the bubble inflated by the jets. 



3 CONCLUSIONS 

Our results provide strong support to the idea that magnetic 
fields can play a crucial role in driving powerful GRB jets 
and associated stellar explosions not only in the magnetar 
model but also in the collapsar model. The main energy 
source for the jets and explosions in our simulations is the 
rotational energy of black hole and it is released via the 
Blandford-Znajek mechanism. The measured rate of energy 
release, E > 2 x lO^^ergs"^, can explain the energetics of 
even the shortest of long duration GRBs. The fact that the 
rotational energy of black hole, E^h — few x lO^^erg, exceeds 
the typical explosion values derived from observations, E ~ 
lO^^erg, suggests a self-regulating process in which the black 
hole activity ceases when the blast wave terminates further 
mass supply to the accretion disk. The full details of the 
simulations together with the results of parameter study will 
be presented elsewhere. 
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